Practica Fundamentos
Ana Moscoso Pérez , Gonzalo Gutiérrez Acosta, Juan Francisco Romero Moreno
17/12/2021
1. Introducción al estudio
Hemos escogido un dataset que recoge ciertas características de unas casas en el área del condado de King, del estado de Whasington junto al valor por el que se vendieron. En el conjunto tenemos un total de 21613 observaciones, es decir viviendas.
Nuestro objetivo será usar esta muestra para poder predecir el valor de una casa en este área en funcion de las variables estudiadas.
Las variables que nos ocupan y que se presentan en el dataset son:
- id – Identificador único de cada casa
- date – Fecha en la que se vendió la casa
- price – Precio
- bedrooms – Número de habitaciones
- bathrooms – Número de baños
- sqft_living – Pies cuadrados útiles de la casa
- sqft_lot – Pies cuadrados del terreno en el que se sitúa la casa
- floors – Número de plantas
- waterfront – Si el apartamento tiene o no vistas al mar
- view – Cómo de buena son las vistas desde la propiedad
- condition – Condición de la casa
- grade – Nivel de construcción y diseño de la casa
- sqft_above – Pies cuadrados del interior de la casa sin contar el sótano
- sqft_basement – Pies cuadrados del sótano
- yr_built – Año en que fue construida la casa inicialmente
- yr_renovated – Año en el que la casa fue renovada
- zipcode – Cógido referente a la zona donde se encuentra la casa
- lat - Latitud
- long - Longitud
- sqft_living15 – media de los pies cuadrados útiles de las 15 casas más cercanas
- sqft_lot15 – media de los pies cuadrados del terreno de las 15 casas más cercanas
2. Tranformaciones en el tipo de variables
str(datos)## 'data.frame': 21613 obs. of 21 variables:
## $ id : num 7.13e+09 6.41e+09 5.63e+09 2.49e+09 1.95e+09 ...
## $ date : chr "20141013T000000" "20141209T000000" "20150225T000000" "20141209T000000" ...
## $ price : num 221900 538000 180000 604000 510000 ...
## $ bedrooms : int 3 3 2 4 3 4 3 3 3 3 ...
## $ bathrooms : num 1 2.25 NA 3 2 4.5 NA 1.5 NA 2.5 ...
## $ sqft_living : int 1180 2570 770 1960 1680 5420 1715 1060 1780 1890 ...
## $ sqft_lot : int 5650 7242 10000 5000 8080 101930 6819 9711 7470 6560 ...
## $ floors : num 1 2 NA 1 1 1 NA 1 NA 2 ...
## $ waterfront : int 0 0 0 0 0 0 0 0 0 0 ...
## $ view : int 0 0 NA 0 0 0 NA 0 NA 0 ...
## $ condition : int 3 3 3 5 3 3 3 3 3 3 ...
## $ grade : int 7 7 6 7 8 11 7 7 7 7 ...
## $ sqft_above : int 1180 2170 770 1050 1680 3890 1715 1060 1050 1890 ...
## $ sqft_basement: int 0 400 0 910 0 1530 0 0 730 0 ...
## $ yr_built : int 1955 1951 1933 1965 1987 2001 1995 1963 1960 2003 ...
## $ yr_renovated : int 0 1991 0 0 0 0 0 0 0 0 ...
## $ zipcode : int 98178 98125 NA 98136 98074 98053 NA 98198 NA 98038 ...
## $ lat : num 47.5 47.7 NA 47.5 47.6 ...
## $ long : num -122 -122 NA -122 -122 ...
## $ sqft_living15: int 1340 1690 2720 1360 1800 4760 2238 1650 1780 2390 ...
## $ sqft_lot15 : int 5650 7639 8062 5000 7503 101930 6819 9711 8113 7570 ...Podemos observar que hay variables que no están codificadas correctamente:
- La variable “id” aparece como numérica,“num”, cuando realmente debería ser de tipo carácter,“chr”.
- La variable “date” aparece como “chr”, cuando debería ser “Date”.
- Las variables “waterfront” y “zipcode” deberían ser de tipo “factor” en vez de enteros “int”.
Por lo tanto, para modificar los tipos de variables hacemos lo siguiente:
datos$id=as.character(datos$id)
datos$date = as.Date(datos$date, "%Y%m%dT%H%M%S")
#datee=datos$date #lo usaremos más adelante, en el gráfico para estudiar la variable date
DATOS<- read.csv("kc_house_data.csv")
datee= as.Date(DATOS$date, "%Y%m%dT%H%M%S")
datos$waterfront = as.factor(datos$waterfront)
datos$zipcode = as.factor(datos$zipcode)
datos$condition=as.factor(datos$condition)
datos$date=as.factor(datos$date)
datos$view=as.factor(datos$view)
str(datos)## 'data.frame': 21613 obs. of 21 variables:
## $ id : chr "7129300520" "6414100192" "5631500400" "2487200875" ...
## $ date : Factor w/ 372 levels "2014-05-02","2014-05-03",..: 165 221 291 221 284 11 57 252 340 306 ...
## $ price : num 221900 538000 180000 604000 510000 ...
## $ bedrooms : int 3 3 2 4 3 4 3 3 3 3 ...
## $ bathrooms : num 1 2.25 NA 3 2 4.5 NA 1.5 NA 2.5 ...
## $ sqft_living : int 1180 2570 770 1960 1680 5420 1715 1060 1780 1890 ...
## $ sqft_lot : int 5650 7242 10000 5000 8080 101930 6819 9711 7470 6560 ...
## $ floors : num 1 2 NA 1 1 1 NA 1 NA 2 ...
## $ waterfront : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ view : Factor w/ 5 levels "0","1","2","3",..: 1 1 NA 1 1 1 NA 1 NA 1 ...
## $ condition : Factor w/ 5 levels "1","2","3","4",..: 3 3 3 5 3 3 3 3 3 3 ...
## $ grade : int 7 7 6 7 8 11 7 7 7 7 ...
## $ sqft_above : int 1180 2170 770 1050 1680 3890 1715 1060 1050 1890 ...
## $ sqft_basement: int 0 400 0 910 0 1530 0 0 730 0 ...
## $ yr_built : int 1955 1951 1933 1965 1987 2001 1995 1963 1960 2003 ...
## $ yr_renovated : int 0 1991 0 0 0 0 0 0 0 0 ...
## $ zipcode : Factor w/ 70 levels "98001","98002",..: 67 56 NA 59 38 30 NA 69 NA 24 ...
## $ lat : num 47.5 47.7 NA 47.5 47.6 ...
## $ long : num -122 -122 NA -122 -122 ...
## $ sqft_living15: int 1340 1690 2720 1360 1800 4760 2238 1650 1780 2390 ...
## $ sqft_lot15 : int 5650 7639 8062 5000 7503 101930 6819 9711 8113 7570 ...Con “str(datos)” podemos ver cómo ya si tenemos las variables bien codificadas.
3. División de los datos es training y testing
Para poder crear nuestro modelo y ver su capacidad de predicción dividiremos el conjunto de datos en dos con un porcentaje de 70% y 30% para el grupo de training y test, respectivamente.
set.seed(12345)
inTraining <- createDataPartition(datos$price,
p = .7, list = FALSE, times = 1)
d_training <- slice(datos, inTraining)
d_testing <- slice(datos, -inTraining)
ntrain=nrow(d_training)4. Análisis exploratorio de datos faltantes: VIM
Veamos cuántas observaciones con datos faltantes tenemos en nuestro dataset:
ntrain-sum(complete.cases(d_training))## [1] 68Obtenemos que hay 68 observaciones de un total de 15131, que tienen datos faltantes. Esto se corresponde con una proporcion de 0.0043 de casos con algún dato faltante. Lo cual al ser menor de 0.03 no tenemos un problema grave. Entre las posibles técnicas a aplicar para la imputación de estos datos faltantes podría ser el asignar la mediana del resto de sus valores conocidos en el caso de variables cualitativas y en el caso de las variables cualitativas asignar la categoría más frecuente de sus conocidas.
Otra opción, debido al bajo número de observaciones con datos faltantes, sería descartar directamente los casos con NA.
Realizaremos unos gráficos para analizar posibles relaciones y estructuras de los valores faltantes.
aggr_plot=aggr(d_training,,col=c("green","red"),numbers=TRUE, sortVars=TRUE, labels=names(d_training),cex.axis=0.7,gap=3,ylab=c("Hist","Pat")) 
##
## Variables sorted by number of missings:
## Variable Count
## view 0.0040314586
## zipcode 0.0022470425
## lat 0.0022470425
## grade 0.0011235212
## waterfront 0.0009913423
## sqft_above 0.0007930738
## id 0.0006608949
## date 0.0006608949
## sqft_living 0.0003304474
## sqft_living15 0.0003304474
## bathrooms 0.0001982685
## floors 0.0001982685
## long 0.0001982685
## price 0.0000000000
## bedrooms 0.0000000000
## sqft_lot 0.0000000000
## condition 0.0000000000
## sqft_basement 0.0000000000
## yr_built 0.0000000000
## yr_renovated 0.0000000000
## sqft_lot15 0.0000000000Los gráficos anteriores nos permiten descubrir rápidamente qué variables de nuestro dataset tienen mayor cantidad de datos faltantes (gráfico a la izquierda) y si puede existir algún patrón de co-ocurrencia en los datos faltantes de varias variables (gráfico a la derecha).
summary(aggr_plot)##
## Missings per variable:
## Variable Count
## id 10
## date 10
## price 0
## bedrooms 0
## bathrooms 3
## sqft_living 5
## sqft_lot 0
## floors 3
## waterfront 15
## view 61
## condition 0
## grade 17
## sqft_above 12
## sqft_basement 0
## yr_built 0
## yr_renovated 0
## zipcode 34
## lat 34
## long 3
## sqft_living15 5
## sqft_lot15 0
##
## Missings in combinations of variables:
## Combinations Count Percent
## 0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0 15063 99.55059150
## 0:0:0:0:0:0:0:0:0:0:0:0:1:0:0:0:0:0:0:0:0 2 0.01321790
## 0:0:0:0:0:0:0:0:0:1:0:0:0:0:0:0:1:1:0:0:0 31 0.20487740
## 0:0:0:0:0:0:0:0:0:1:0:1:0:0:0:0:0:0:0:0:0 17 0.11235212
## 0:0:0:0:0:1:0:0:1:0:0:0:0:0:0:0:0:0:0:1:0 5 0.03304474
## 0:0:0:0:1:0:0:1:0:1:0:0:0:0:0:0:1:1:1:0:0 3 0.01982685
## 1:1:0:0:0:0:0:0:1:1:0:0:1:0:0:0:0:0:0:0:0 10 0.06608949Aquí podemos observar el porcentaje de valores perdidos para cada variable. Dichos porcentajes nos permiten ver que no existe ningun patrón de cocurrencia en los datos faltantes de las variables. Esto es, que no tenemos valores de variables que tienden a suceder conjuntamente.
Debido al bajo número de observaciones con datos faltantes y a que no existe ningún patrón de cocurrencia en los datos faltantes procedemos a eliminarlos de nuestro estudio.
d_training=d_training[complete.cases(d_training)==TRUE,] #eliminamos los datos faltantes
dim(d_training)## [1] 15063 21Ahora trabajaremos con un conjunto de datos con 15063 observaciones.
5. Estudio de las variables y sus correspondientes transformaciones
Price
El histograma es un gráfico que representa de forma bastante precisa la distribución del conjunto de datos pudiéndose observar su dispersión y/o asimetría. Por ello, realizaremos un histograma del precio para ver cómose comporta esta variable.
d_training %>% ggplot( aes(x=price)) + geom_histogram(color="black", fill="yellow", lwd=0.3)+
ggtitle("Histograma del precio") +ylab("Frecuencia")
Observamos que la distribución de la variable precio presenta un alto sesgo positivo, por lo que le vamos a aplicar una transformación logarítmica para ver si así consegimos una distribución más simétrica.
p1=d_training %>% ggplot( aes(x= price)) + geom_histogram(aes(y=..density..), bins=30, colour="black", fill="yellow") +
geom_density(alpha=.5, fill="blue") + ggtitle("Precio")
p2=d_training %>% mutate(log10_p=log(price)) %>% ggplot( aes(x=log10_p)) + geom_histogram(aes(y=..density..), bins=30, colour="black", fill="yellow") +
geom_density(alpha=.5, fill="blue") + ggtitle("Logaritmo del precio")
grid.arrange(p1,p2,nrow=2)
Ahora obtenemos una distribución que a primera vista es como una campana de gauss, de modo que podríamos afirmar que dicha variable sigue una distribución normal. Por ello, a partir de ahora trabajaremos con la variable log(price) en vez de con price.
d_training$price=log(d_training$price)Bedrooms
d1=d_training %>% ggplot( aes(x=bedrooms)) + geom_histogram(colour="black", bins =30,fill="tomato")+
ylab("Frecuencias")
d2=d_training %>% ggplot( aes(x=as.factor(bedrooms), y=price, fill=as.factor(bedrooms))) + geom_boxplot()+
labs(x="bedrooms")+theme(legend.position="none")
d3=ggplot(d_training, aes(x = "", y = bedrooms)) + stat_boxplot(geom = "errorbar", width = 0.3) +
geom_boxplot(fill = "orange", outlier.colour = "red", alpha = 0.9) +
xlab("")
grid.arrange(d1,d2,d3,nrow=2) 
En los gráficos representados podemos observar que según aumenta el número de habitaciones aumenta el precio de las casas. También podemos ver posibles outliers, como las casas con 0 y 6 o más habitaciones. Estudiemoslas.
sum(d_training$bedrooms == 0)## [1] 13Vemos que hay 13 casas con cero habitaciones. Aunque no es frecuente hay casas, llamadas estudios, en las cuales todo está en una misma estancia. Por ello no consideramos estas observaciones como datos erróneos y las seguimos manteniendo en nuestro dataset.
sum(d_training$bedrooms >=6) ## [1] 230Hay 228 casas con 6 y más habitaciones. Al ser un número elevado podemos intuir que tener alto número de habitaciones es algo normal por la zona.
sum(d_training$bedrooms > 8)## [1] 8Sin embargo, solo hay 7 casas con 9 o más habitaciones. Esto ya es algo más peculiar, de modo que necesitarán un estudio más profundo.
En principio comprobamos en el mapa la hubicación de estas casas, por si se diese el caso de que están cerca. Si esto fuese así podríamos pensar que es propio de la zona.
Observamos que las casas no están en un misma zona. Lo siguiente que haremos será ver el número de habitaciones que tienen las casas de alrededor.
Comprobando el número de habitaciones de las casas cercanas a nuestros posibles outliers observamos que estas tienen entre 2 y 3 habitaciones. Por ello, el que haya casas por esas zonas con altos números de dormitorios es un poco extraño. Ya que el número de observaciones que tienen más de 8 habitaciones respecto al total de los datos de entrenamiento es muy poco representativo las eliminaremos del estudio.
d_training = d_training[d_training$bedrooms < 8,]Bathrooms
Observemos los posibles valores que toma esta variable:
table(d_training$bathrooms)##
## 0 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3 3.25 3.5 3.75 4
## 10 4 47 2667 5 1007 2135 1341 1417 3779 807 519 410 514 111 90
## 4.25 4.5 4.75 5 5.25 5.5 5.75 6 6.25 6.5 6.75 7.75
## 56 69 19 14 9 5 3 3 1 2 1 1La variable Bathrooms puede tomar valores decimales de 0.25 en 0.25. Y va desde 0 a 8. Esto es debido a que el número de baños se contabiliza por las piezas y cada baño completo tiene 4 piezas.
- Baño (4 piezas) -> Inodoro, lavabo, bañera y ducha.(1 unidad) - Baño completo
- Baño (3 piezas) -> Inodoro, lavabo y ducha. (0.75 unidad) - Baño con ducha
- Baño (2 piezas) -> Inodoro y lavabo. (0.5 unidad) - Aseo
- Baño (1 pieza) -> Inodoro (0.25 unidad) - Aseo
Vamos a agrupar los baños del siguiente modo:
- Valores: 0.75, 1 == 1 baño
- Valores: 1.25, 1.5, 1.75, 2 == 2 baños
- Valores: 2.25, 2.5, 2.75, 3 == 3 baños
Y así hasta tener 8 baños.
Sin embargo, las casas cuya variable baño tome un valor menor de 0.75 no las consideraremos, ya que una casa que no tenga como mínimo un ducha o un lavabo no es algo lógico. Dependiendo del número de casos en los que ocurra esto, los tomaremos como datos faltantes o simplemente errores de mediciones y si no son muchas las elimiaremos del estudio.
sum(d_training$bathrooms<0.75)## [1] 14Obtenemos 14 casas. Que es una número irrelevante frente al total de casos que tenemos. Por lo que eliminaremos estos casos del estudio.
d_training=d_training[-which(d_training$bathrooms<0.75),]A continuación lo que haremos será cambiar la variable para que considere número de baños en vez de estancias.
d_training$bathrooms=ceiling(d_training$bathrooms)ggplot(d_training, aes(x = "", y = bathrooms)) +
stat_boxplot(geom = "errorbar", width = 0.3) +
geom_boxplot(fill = "orange", outlier.colour = "red", alpha = 0.9) +
ggtitle("Boxplot bathrooms") +
xlab("") +
coord_flip() 
table(d_training$bathrooms) ##
## 1 2 3 4 5 6 7 8
## 2714 4488 6522 1125 158 20 4 1Podemos observar que hay dos casas con 8 baños y 4 casas con 7 baños. Debido a que esto no es algo común las analizaremos.´
Estudiamos las casas con 7 baños.
d_training[d_training$bathrooms==7,] ## id date price bedrooms bathrooms sqft_living sqft_lot
## 5667 1924059029 2014-06-17 15.35624 5 7 9640 13068
## 10182 2303900035 2014-06-11 14.87607 5 7 8670 64033
## 14398 424069279 2015-03-28 13.98102 6 7 6260 10955
## 15049 2524069097 2014-05-09 14.62149 5 7 7270 130017
## floors waterfront view condition grade sqft_above sqft_basement yr_built
## 5667 1 1 4 3 12 4820 4820 1983
## 10182 2 0 4 3 13 6120 2550 1965
## 14398 2 0 0 3 11 4840 1420 2007
## 15049 2 0 0 3 12 6420 850 2010
## yr_renovated zipcode lat long sqft_living15 sqft_lot15
## 5667 2009 98040 47.5570 -122.210 3270 10454
## 10182 2003 98177 47.7295 -122.372 4140 81021
## 14398 0 98075 47.5947 -122.039 2710 12550
## 15049 0 98027 47.5371 -121.982 1800 44890Nos concuerda el número de baños con los pies cuadrados de la casa, por lo que las conservaremos en nuestro conjunto de datos.
Ahora estudiaremos la casa con 8 baños:
d_training[d_training$bathrooms==8,]## id date price bedrooms bathrooms sqft_living sqft_lot
## 6467 9208900037 2014-09-19 15.74486 6 8 9890 31374
## floors waterfront view condition grade sqft_above sqft_basement yr_built
## 6467 2 0 4 3 13 8860 1030 2001
## yr_renovated zipcode lat long sqft_living15 sqft_lot15
## 6467 0 98039 47.6305 -122.24 4540 42730Observamos un caso que tiene más baños que habitaciones y otro caso que tiene demasiados baños y habitaciones respecto a los pies cuadrados de la casa. Lo cual no tiene mucho sentido. Por lo tanto, eliminamos estas dos observaciones del estudio.
d_training=d_training[-which(d_training$bathrooms==8),]Agrupamos las casas según el número de baños y realizamos un boxplot para ver si hay diferencias o no en el precio.
d_training %>% ggplot(aes(x=as.factor(d_training$bathrooms), y=price, fill=as.factor(d_training$bathrooms))) + geom_boxplot()+
labs(x="bathrooms_group")+theme(legend.position="none")
Se puede ver claramente como según van aumentando el número de baños, aumenta el precio de la casa.
Sqft_living
Veamos su comportamiento frente al precio:
l1<-d_training %>% ggplot(aes(x=sqft_living)) +
geom_histogram(aes(y=..density..), bins=30, colour="black", fill="yellow") +
geom_density(alpha=.3, fill="blue")
l2<-d_training %>% ggplot(aes(sqft_living, price)) +
geom_point(alpha = 0.5,col="blue") +
geom_smooth(se = F, method = "lm", color = "red") +
scale_y_continuous(breaks = seq(0,8000000, by = 1000000))
grid.arrange(l1,l2, nrow=1)
Viendo el scaterplot de los pies cuadrados de la casa frente al precio (gráfico derecho) podemos ver que hay una relación creciente, a medida que aumentan los pies cuadrados de las casas aumenta el precio. Además, en el gráfico izquierdo podemos observar como la distribución de la variable presenta un alto sesgo positivo. Debido a que sería conveniente transformarla para que la distribución de valores fuese más homogénea, le aplicamos una transformación logarítmica.
d_training$sqft_living= log(d_training$sqft_living)
l1log <- ggplot(d_training, aes(x=sqft_living)) + geom_histogram(aes(y=..density..), bins=30, colour="black", fill="white") + geom_density(alpha=.3, fill="#E1AF00")
l2log <- ggplot(d_training, aes(sqft_living, price)) +
geom_point(alpha = 0.5,col="blue") +
geom_smooth(se = F, method = "lm", color = "red") +
scale_y_continuous(breaks = seq(0,8000000, by = 1000000))
grid.arrange(l1log,l2log, nrow=1)
Se puede observar cómo ahora los datos se han normalizado (la distribución presenta la forma de una campana de Gauss).
Basement y above
A continuación vamos a comprobar que la variable above más la variable basement coincide con la variable living. A partir de aquí, vamos a deducir que casas tienen sótano y si el tenerlo o no afecta en el precio.
Total=log(d_training$sqft_above+d_training$sqft_basement)
all(Total==d_training$sqft_living) ## [1] TRUE#comprobamos que efectivamente sqft_above+sqft_basement=sqft_living
a=cbind(Total,living=d_training$sqft_living,liv15=log(d_training$sqft_living15),lot=log(d_training$sqft_lot),lot15=(d_training$sqft_lot15))
head(a)## Total living liv15 lot lot15
## [1,] 7.073270 7.073270 7.200425 8.639411 5650
## [2,] 7.851661 7.851661 7.432484 8.887653 7639
## [3,] 7.580700 7.580700 7.215240 8.517193 5000
## [4,] 7.426549 7.426549 7.495542 8.997147 7503
## [5,] 8.597851 8.597851 8.468003 11.532042 101930
## [6,] 7.544332 7.544332 7.779049 8.788746 7570d_training$basement=ifelse(d_training$sqft_living-log(d_training$sqft_above)==0,"0","1") #0 es que no tiene sotano
d_training$basement=as.factor(d_training$basement)
ggplot(d_training, aes(x=basement, y=price, fill=basement)) + geom_boxplot()+
labs(x="Tener sótano (1) o no tenerlo (0)")+
scale_y_continuous(labels = scales::dollar,n.breaks = 15)
Se ve un pequeño incremento en el precio en las casas que tienen sótano. Pero no es perceptible, por lo que no consideraremos el tener o no sótano relevante en el modelo.
Living15
Para comprobar si podemos explicar la variable sqft_living15 a partir de sqft_living realizaremos una recta de regresión:
head(a)## Total living liv15 lot lot15
## [1,] 7.073270 7.073270 7.200425 8.639411 5650
## [2,] 7.851661 7.851661 7.432484 8.887653 7639
## [3,] 7.580700 7.580700 7.215240 8.517193 5000
## [4,] 7.426549 7.426549 7.495542 8.997147 7503
## [5,] 8.597851 8.597851 8.468003 11.532042 101930
## [6,] 7.544332 7.544332 7.779049 8.788746 7570cbind(mean(a[,3]-a[,2]),mean(a[,5]-a[,4]))## [,1] [,2]
## [1,] -0.009071659 12714.12summary(lm(log(d_training$sqft_living15)~d_training$sqft_living, data=d_training))##
## Call:
## lm(formula = log(d_training$sqft_living15) ~ d_training$sqft_living,
## data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.65892 -0.13088 0.00365 0.14332 1.13938
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.211374 0.031654 101.5 <2e-16 ***
## d_training$sqft_living 0.573466 0.004186 137.0 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2183 on 15029 degrees of freedom
## Multiple R-squared: 0.5553, Adjusted R-squared: 0.5553
## F-statistic: 1.877e+04 on 1 and 15029 DF, p-value: < 2.2e-16Realizamos un contraste de hipótesis en el cual la hipótesis nula es que la variable no es relevante para el modelo. Esto es, no podemos explicar la variable respuesta (sqft_living15) en función de la explicativa (sqft_living). Al obtener un p.valor menor que 0.05 rechazamos H0 y por consiguiente podemos explicar living15 a partir de living. (Aún así no obtenemos un R^2 ajustado bueno que digamos).
Veamos que variable de entre estas dos explica mejor el precio:
summary(lm(d_training$price~d_training$sqft_living15, data=d_training)) ##
## Call:
## lm(formula = d_training$price ~ d_training$sqft_living15, data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.88918 -0.30016 -0.00795 0.26023 1.81985
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.211e+01 1.037e-02 1167.04 <2e-16 ***
## d_training$sqft_living15 4.736e-04 4.928e-06 96.11 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4145 on 15029 degrees of freedom
## Multiple R-squared: 0.3806, Adjusted R-squared: 0.3806
## F-statistic: 9237 on 1 and 15029 DF, p-value: < 2.2e-16summary(lm(d_training$price~d_training$sqft_living, data=d_training)) ##
## Call:
## lm(formula = d_training$price ~ d_training$sqft_living, data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10747 -0.29387 0.01419 0.25964 1.31752
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.749564 0.056446 119.6 <2e-16 ***
## d_training$sqft_living 0.834292 0.007464 111.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3892 on 15029 degrees of freedom
## Multiple R-squared: 0.4539, Adjusted R-squared: 0.4539
## F-statistic: 1.249e+04 on 1 and 15029 DF, p-value: < 2.2e-16Prediciendo el precio según living obtenemos un R² un poco mayor. Luego entre living y living15 nos quedaremos con living.
Sqft_lot
sl0=ggplot(d_training, aes(x = "", y = sqft_lot)) +
stat_boxplot(geom = "errorbar", width = 0.3) +
geom_boxplot(fill = "orange", outlier.colour = "red",alpha = 0.9) +
xlab("") + coord_flip()
sl1=ggplot(d_training, aes(x=sqft_lot)) + geom_histogram(aes(y=..density..), bins=30, colour="black", fill="tomato")
sl2=ggplot(d_training, aes(sqft_lot, price)) +
geom_point(alpha = 0.5) +
geom_smooth(se = F, method = "lm", color = "blue") +
scale_y_continuous(breaks = seq(0,8000000, by = 1000000))
grid.arrange(sl0,sl1,sl2, nrow=2)
Esta variable no parece que tenga una relación lineal muy clara con price, por lo que no vemos necesario hacer ninguna transformación ya que no la usaremos para la implementación del modelo. Haciendo el estudio de manera análoga para lot15 llegamos a la misma conclusión.
Floors
table(d_training$floors)##
## 1 1.5 2 2.5 3 3.5
## 7397 1328 5755 115 434 2No se puede tener un número de pisos que no sea entero. Debido a que no tenemos muy claro el significado de estos decimales, supondremos que indica que la planta está inacabada o que hay espacio en la misma que no se está usando. Atendiendo a esta consideración sobre el número de plantas truncaremos la variable.
d_training$floors=trunc(d_training$floors)De esta forma nos quedamos con 3 posibles pisos. Veamos si esto influye o no en el precio de las casas.
f1=ggplot(d_training, aes(x=as.factor(floors), y=price, fill=as.factor(floors))) + geom_boxplot()+
labs(x="Waterfront")+
scale_y_continuous(labels = scales::dollar,n.breaks = 15)
f2=ggplot(d_training, aes(x=as.factor(floors))) + geom_histogram(colour="black", stat ="count",fill="yellow")+
labs(x="Condición",y="Frecuencia")
grid.arrange(f2,f1, nrow=1)
Vemos que el numero de plantas no influye en el precio de las casas.
Waterfront
La codificación respecto si tiene o no vistas al mar es 0 y 1, respectivamente.
ggplot(d_training, aes(x=waterfront, y=price, fill=waterfront)) + geom_boxplot()+
labs(x="Waterfront")+
scale_y_continuous(labels = scales::dollar,n.breaks = 15)
Observamos un incremento notorio en el precio de las casas que tienen vistas al mar.
Yr_renovated
Vamos a hacer un cambio a esta variable. Y la cambiaremos a sí o no en función de si ha sido o no renovada. Codificamos 0 como no renovada y 1 como renovada, y comprobamos si el que la casa haya sido o no renovada influye en el precio.
d_training$yr_renovated=ifelse(d_training$yr_renovated==0,"0", "1")
d_training$yr_renovated=as.factor(d_training$yr_renovated)
ggplot(d_training, aes(x=yr_renovated, y=price, fill=yr_renovated)) + geom_boxplot()+
labs(x="No renovada VS Renovada")+theme(legend.position="none")+
scale_y_continuous(labels = scales::dollar,n.breaks = 15)
Veamos que sí influye en el precio el que la casa haya sido o no renovada.
Yr_built
Las casas se construyeron entre 1900 y 2015. Veamos si influye esto en el precio
summarise(d_training, añoMin=min(d_training$yr_built),añoMax= max(d_training$yr_built))## añoMin añoMax
## 1 1900 2015ggplot(d_training, aes(x=yr_built)) + geom_histogram(colour="black",bins=30, fill="#E1AF00")+
labs(title="Histograma del precio y año de construcción", x="Año construida",y="Precio")
Vemos un claro aumento en el precio de las casas construidas en los últimos años.
Condition
table(d_training$condition)##
## 1 2 3 4 5
## 23 119 9767 3951 1171c1=ggplot(d_training, aes(x=condition, y=price, fill=condition)) + geom_boxplot()+
ylab("Precio")+theme(legend.position="none")+
scale_y_continuous(labels = scales::dollar,n.breaks = 15)
c2=ggplot(d_training, aes(x=condition)) + geom_histogram(colour="black", stat ="count",fill="yellow")+
ylab("Frecuencia")
grid.arrange(c1,c2, nrow=1)
Observamos que sí es una variable significativa.
View
v1=ggplot(d_training, aes(x=view, y=price, fill=view)) + geom_boxplot()+
ylab("Precio")+
scale_y_continuous(labels = scales::dollar,n.breaks = 15)
v2=ggplot(d_training, aes(x=view)) + geom_histogram(colour="black", stat ="count",fill="yellow")+
ylab("Frecuencia")
grid.arrange(v1,v2, nrow=1)
Observamos que si la casa tiene mejores vistas su precio es mayor. Sin embargo, no hay muchas casas con buenas vistas.
Grade
La variable grade mide el nivel de construcción y diseño de la casa en un rango del 1 al 13. Entre los valores 1 y 3 indica que no llega a la construcción y el diseño de edificios, 7 que tiene un nivel promedio, y entre 11 y 13 indica que tiene un alto nivel de calidad de construcción y diseño.
En función de la descripción de la variable tomada, vamos a agruparla en 4 para reducir el número de sus posibles niveles.
- Inacabada = 0 a 3
- Aceptable = 4 a 7
- Buena = 8 a 10
- Excelente = 11 a 13
d_training$grade <- cut(d_training$grade, breaks = c(0,3,7,10,13), labels = c("Inacabada","Aceptable","Buena","Excelente"))
ggplot(d_training, aes(x=grade, y=price, fill=grade)) + geom_boxplot()+
labs(x="Grado de construcción",y="Precio")+
scale_y_continuous(labels = scales::dollar,n.breaks = 15)
El grado de construcción sí que afecta al precio, y bastante como podemos observar.
Zipcode
Realizamos un boxplot diferenciando con colores los distintos valores de zipcode:
ggplot(d_training, aes(x=zipcode, y=price, fill=zipcode)) + geom_boxplot()+
ylab("Precio")+ guides(fill=FALSE) +
scale_y_continuous(labels = scales::dollar,n.breaks = 15) +
theme(
axis.text.x = element_text(angle=90,size=7)
)
A simple vista no se observa ningún patrón de precios dependiendo del área donde esté ubicada la casa.
6. Ajuste, interpretación y diagnosis del modelo de regresión lineal
Nuestro objetivo es predecir el precio de las casas en el condado de King y comprender qué factores son los responsables de un valor de propiedad más alto.
Antes de nada eliminaremos del conjunto de dato la variable “identificador” ya que no es relevante para predecir el precio de la vivienda. No es más que un identificador de la casa, como sus nombre por así decirlo .
d_training=select(d_training,-id)A continuación ajustaremos una recta de regresión para predecir la variable respuesta precio en función de las variables explicativas del modelo.
Veamos en principio la colinealidad entre las varibles cuantitativas
d_train_cuant=select(d_training,-date,-waterfront,-view,-condition,-basement,-grade,-yr_built,-yr_renovated,-zipcode) #seleccionamos las variables cuantitativas
correlaciones= cor(d_train_cuant) #correlación entre las variables cuantitativas
correlaciones## price bedrooms bathrooms sqft_living sqft_lot
## price 1.00000000 0.355195293 0.52864672 0.67374379 0.090451314
## bedrooms 0.35519529 1.000000000 0.50620174 0.64462664 0.027010099
## bathrooms 0.52864672 0.506201745 1.00000000 0.74934905 0.067675948
## sqft_living 0.67374379 0.644626641 0.74934905 1.00000000 0.134933646
## sqft_lot 0.09045131 0.027010099 0.06767595 0.13493365 1.000000000
## floors 0.28801889 0.163688858 0.55381091 0.36347719 -0.009682656
## sqft_above 0.60316447 0.499271684 0.66635231 0.83440895 0.171245781
## sqft_basement 0.31084025 0.301591608 0.25086854 0.41639899 -0.001037189
## lat 0.45036426 -0.008473969 0.01652289 0.03893912 -0.094040091
## long 0.04919025 0.144571652 0.24061825 0.25933751 0.227252958
## sqft_living15 0.61696436 0.410253085 0.56037916 0.73464487 0.135475873
## sqft_lot15 0.08100238 0.026356736 0.06888082 0.15188713 0.696687244
## floors sqft_above sqft_basement lat long
## price 0.288018888 0.603164473 0.310840247 0.450364264 0.04919025
## bedrooms 0.163688858 0.499271684 0.301591608 -0.008473969 0.14457165
## bathrooms 0.553810915 0.666352306 0.250868537 0.016522892 0.24061825
## sqft_living 0.363477193 0.834408954 0.416398988 0.038939116 0.25933751
## sqft_lot -0.009682656 0.171245781 -0.001037189 -0.094040091 0.22725296
## floors 1.000000000 0.518092745 -0.242020536 0.032012208 0.15984579
## sqft_above 0.518092745 1.000000000 -0.062427518 -0.002549824 0.34877063
## sqft_basement -0.242020536 -0.062427518 1.000000000 0.112435799 -0.15086633
## lat 0.032012208 -0.002549824 0.112435799 1.000000000 -0.14519166
## long 0.159845794 0.348770625 -0.150866335 -0.145191661 1.00000000
## sqft_living15 0.291342424 0.728462262 0.197486776 0.045002865 0.33933155
## sqft_lot15 -0.020755199 0.179429883 0.012308039 -0.101471044 0.25137265
## sqft_living15 sqft_lot15
## price 0.61696436 0.08100238
## bedrooms 0.41025309 0.02635674
## bathrooms 0.56037916 0.06888082
## sqft_living 0.73464487 0.15188713
## sqft_lot 0.13547587 0.69668724
## floors 0.29134242 -0.02075520
## sqft_above 0.72846226 0.17942988
## sqft_basement 0.19748678 0.01230804
## lat 0.04500286 -0.10147104
## long 0.33933155 0.25137265
## sqft_living15 1.00000000 0.18380430
## sqft_lot15 0.18380430 1.00000000corrplot(correlaciones,type="upper")
Correlaciones con la variable respuesta:
correlaciones[,1] ## price bedrooms bathrooms sqft_living sqft_lot
## 1.00000000 0.35519529 0.52864672 0.67374379 0.09045131
## floors sqft_above sqft_basement lat long
## 0.28801889 0.60316447 0.31084025 0.45036426 0.04919025
## sqft_living15 sqft_lot15
## 0.61696436 0.08100238Con respecto al precio vemos alta correlación con el sqft_living y sqft_above. De modo que en el modelo futuro se espera que estas variables tengan mayor peso.
det(correlaciones)## [1] 0.0004694057El determinante de la matriz de correlaciones entre las variables explicativas es cero. Esto implica que si consideramos todas estas variables tendremos un problema de multicolinealidad. Por ello, para construir la recta de regresión realizaremos una selección de variables a mano.
Según el estudio visual de las variables más la matriz de correlaciones tomamos en principio las siguiente selección de variables:
model1 <- lm(price~bedrooms+bathrooms+sqft_living+waterfront+view+condition+
grade+yr_renovated+yr_built+lat, data=d_training)
summary(model1)##
## Call:
## lm(formula = price ~ bedrooms + bathrooms + sqft_living + waterfront +
## view + condition + grade + yr_renovated + yr_built + lat,
## data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.260 -0.175 0.000 0.166 1.228
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.525e+01 9.080e-01 -60.854 < 2e-16 ***
## bedrooms -3.877e-02 3.387e-03 -11.448 < 2e-16 ***
## bathrooms 6.125e-02 4.250e-03 14.412 < 2e-16 ***
## sqft_living 5.826e-01 9.878e-03 58.978 < 2e-16 ***
## waterfront1 3.549e-01 3.237e-02 10.963 < 2e-16 ***
## view1 1.976e-01 1.783e-02 11.084 < 2e-16 ***
## view2 1.460e-01 1.093e-02 13.362 < 2e-16 ***
## view3 2.308e-01 1.501e-02 15.373 < 2e-16 ***
## view4 3.040e-01 2.278e-02 13.346 < 2e-16 ***
## condition2 8.360e-02 6.220e-02 1.344 0.178935
## condition3 1.915e-01 5.716e-02 3.350 0.000811 ***
## condition4 2.388e-01 5.715e-02 4.179 2.94e-05 ***
## condition5 2.914e-01 5.752e-02 5.066 4.10e-07 ***
## gradeAceptable -6.251e-01 2.732e-01 -2.288 0.022156 *
## gradeBuena -3.898e-01 2.733e-01 -1.426 0.153851
## gradeExcelente 2.333e-02 2.739e-01 0.085 0.932135
## yr_renovated1 6.725e-02 1.187e-02 5.667 1.48e-08 ***
## yr_built -2.630e-03 1.094e-04 -24.041 < 2e-16 ***
## lat 1.458e+00 1.668e-02 87.397 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2729 on 15012 degrees of freedom
## Multiple R-squared: 0.7319, Adjusted R-squared: 0.7316
## F-statistic: 2277 on 18 and 15012 DF, p-value: < 2.2e-16Obtenemos unos coeficientes lógicos y un R² ajustado de 0.7316. En cambio indica que a mayor número de habitaciones menor precio, lo cual es un poco raro.
En general es buen modelo, aunque la variable grade no parece influir mucho. Vemos que pasa si no la introducimos:
model1b=lm(price~bedrooms+bathrooms+sqft_living+waterfront+view+condition+
yr_renovated+yr_built+lat, data=d_training) #model1 sin el grade
summary(model1b) ##
## Call:
## lm(formula = price ~ bedrooms + bathrooms + sqft_living + waterfront +
## view + condition + yr_renovated + yr_built + lat, data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23325 -0.18387 -0.00279 0.17633 1.25404
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.263e+01 9.153e-01 -68.426 < 2e-16 ***
## bedrooms -6.310e-02 3.571e-03 -17.672 < 2e-16 ***
## bathrooms 9.114e-02 4.483e-03 20.329 < 2e-16 ***
## sqft_living 7.487e-01 9.794e-03 76.440 < 2e-16 ***
## waterfront1 3.585e-01 3.454e-02 10.379 < 2e-16 ***
## view1 2.222e-01 1.902e-02 11.679 < 2e-16 ***
## view2 1.874e-01 1.163e-02 16.117 < 2e-16 ***
## view3 2.874e-01 1.597e-02 17.997 < 2e-16 ***
## view4 3.727e-01 2.425e-02 15.366 < 2e-16 ***
## condition2 6.170e-02 6.640e-02 0.929 0.352731
## condition3 1.766e-01 6.102e-02 2.894 0.003809 **
## condition4 2.147e-01 6.101e-02 3.520 0.000433 ***
## condition5 2.601e-01 6.140e-02 4.235 2.30e-05 ***
## yr_renovated1 6.729e-02 1.265e-02 5.318 1.06e-07 ***
## yr_built -1.907e-03 1.145e-04 -16.658 < 2e-16 ***
## lat 1.547e+00 1.766e-02 87.572 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2913 on 15015 degrees of freedom
## Multiple R-squared: 0.6944, Adjusted R-squared: 0.6941
## F-statistic: 2275 on 15 and 15015 DF, p-value: < 2.2e-16model1c <- lm(price~bathrooms+sqft_living+waterfront+view+condition+
grade+yr_renovated+yr_built+lat, data=d_training) #modelo1 sin bedroom
summary(model1c)##
## Call:
## lm(formula = price ~ bathrooms + sqft_living + waterfront + view +
## condition + grade + yr_renovated + yr_built + lat, data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.26112 -0.17537 0.00039 0.16626 1.23274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.563e+01 9.113e-01 -61.042 < 2e-16 ***
## bathrooms 5.489e-02 4.232e-03 12.971 < 2e-16 ***
## sqft_living 5.272e-01 8.648e-03 60.963 < 2e-16 ***
## waterfront1 3.681e-01 3.249e-02 11.330 < 2e-16 ***
## view1 2.051e-01 1.789e-02 11.463 < 2e-16 ***
## view2 1.529e-01 1.096e-02 13.949 < 2e-16 ***
## view3 2.418e-01 1.505e-02 16.071 < 2e-16 ***
## view4 3.141e-01 2.286e-02 13.742 < 2e-16 ***
## condition2 8.534e-02 6.246e-02 1.366 0.171910
## condition3 1.903e-01 5.741e-02 3.315 0.000919 ***
## condition4 2.364e-01 5.739e-02 4.118 3.84e-05 ***
## condition5 2.887e-01 5.777e-02 4.998 5.86e-07 ***
## gradeAceptable -6.527e-01 2.744e-01 -2.379 0.017370 *
## gradeBuena -4.079e-01 2.745e-01 -1.486 0.137278
## gradeExcelente 1.995e-02 2.751e-01 0.073 0.942200
## yr_renovated1 7.311e-02 1.191e-02 6.140 8.43e-10 ***
## yr_built -2.503e-03 1.093e-04 -22.899 < 2e-16 ***
## lat 1.467e+00 1.673e-02 87.695 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2741 on 15013 degrees of freedom
## Multiple R-squared: 0.7296, Adjusted R-squared: 0.7293
## F-statistic: 2383 on 17 and 15013 DF, p-value: < 2.2e-16El modelo llamado “model1b” indica que mientras más habitaciones más baratas las casas… sospechoso y un R² = 0.6941. EL model1c presenta en general unos coeficientes lógicos, aunque los de condi2 y vieew1 nos chocan un poco y un R²=0.7293.
model1cc <- lm(price~bathrooms+sqft_living+waterfront+view+condition+
yr_renovated+yr_built+lat, data=d_training) # modelo1c sin grade
summary(model1cc)##
## Call:
## lm(formula = price ~ bathrooms + sqft_living + waterfront + view +
## condition + yr_renovated + yr_built + lat, data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23353 -0.18880 -0.00439 0.18000 1.26612
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.376e+01 9.225e-01 -69.111 < 2e-16 ***
## bathrooms 8.249e-02 4.503e-03 18.320 < 2e-16 ***
## sqft_living 6.671e-01 8.727e-03 76.438 < 2e-16 ***
## waterfront1 3.806e-01 3.487e-02 10.915 < 2e-16 ***
## view1 2.363e-01 1.920e-02 12.306 < 2e-16 ***
## view2 2.015e-01 1.172e-02 17.195 < 2e-16 ***
## view3 3.094e-01 1.608e-02 19.236 < 2e-16 ***
## view4 3.940e-01 2.447e-02 16.098 < 2e-16 ***
## condition2 6.313e-02 6.708e-02 0.941 0.346670
## condition3 1.737e-01 6.165e-02 2.817 0.004847 **
## condition4 2.091e-01 6.164e-02 3.392 0.000695 ***
## condition5 2.536e-01 6.204e-02 4.087 4.39e-05 ***
## yr_renovated1 7.719e-02 1.277e-02 6.045 1.53e-09 ***
## yr_built -1.645e-03 1.147e-04 -14.344 < 2e-16 ***
## lat 1.568e+00 1.780e-02 88.111 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2943 on 15016 degrees of freedom
## Multiple R-squared: 0.6881, Adjusted R-squared: 0.6878
## F-statistic: 2366 on 14 and 15016 DF, p-value: < 2.2e-16Obtenemos un R² de 0.6878, lo cual empeoramos. Rechazamos esta selección de variables.
model1a=(lm(price~sqft_living+waterfront+view+condition+
+grade+yr_renovated+yr_built+lat, data=d_training)) #model1 quitando bathrooms y bedrooms
summary(model1a)#coeficientes lógicos, aunque cond2 y view1 no del todo y R²= 0.7263##
## Call:
## lm(formula = price ~ sqft_living + waterfront + view + condition +
## +grade + yr_renovated + yr_built + lat, data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.20412 -0.17392 0.00029 0.16770 1.22714
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.767e+01 9.026e-01 -63.896 < 2e-16 ***
## sqft_living 5.926e-01 7.061e-03 83.926 < 2e-16 ***
## waterfront1 3.720e-01 3.267e-02 11.386 < 2e-16 ***
## view1 2.065e-01 1.799e-02 11.478 < 2e-16 ***
## view2 1.555e-01 1.102e-02 14.111 < 2e-16 ***
## view3 2.467e-01 1.513e-02 16.312 < 2e-16 ***
## view4 3.137e-01 2.299e-02 13.646 < 2e-16 ***
## condition2 8.442e-02 6.281e-02 1.344 0.17898
## condition3 1.890e-01 5.772e-02 3.275 0.00106 **
## condition4 2.348e-01 5.771e-02 4.068 4.77e-05 ***
## condition5 2.961e-01 5.809e-02 5.098 3.48e-07 ***
## gradeAceptable -6.899e-01 2.759e-01 -2.501 0.01240 *
## gradeBuena -4.360e-01 2.760e-01 -1.580 0.11419
## gradeExcelente 6.933e-03 2.766e-01 0.025 0.98001
## yr_renovated1 9.004e-02 1.190e-02 7.566 4.07e-14 ***
## yr_built -1.922e-03 1.002e-04 -19.170 < 2e-16 ***
## lat 1.479e+00 1.680e-02 88.057 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2756 on 15014 degrees of freedom
## Multiple R-squared: 0.7266, Adjusted R-squared: 0.7263
## F-statistic: 2493 on 16 and 15014 DF, p-value: < 2.2e-16model2 =lm(price~bathrooms+sqft_living+waterfront+view+condition+
grade+yr_renovated+lat, data=d_training) #modelo1 sin bedroom ni yr_built
summary(model2)#los coeficientes del grade no son lógicos##
## Call:
## lm(formula = price ~ bathrooms + sqft_living + waterfront + view +
## condition + grade + yr_renovated + lat, data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.22004 -0.17945 -0.00649 0.16520 1.28628
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -64.813514 0.832423 -77.861 < 2e-16 ***
## bathrooms 0.015177 0.003927 3.865 0.000111 ***
## sqft_living 0.546092 0.008757 62.362 < 2e-16 ***
## waterfront1 0.368346 0.033051 11.145 < 2e-16 ***
## view1 0.227822 0.018172 12.537 < 2e-16 ***
## view2 0.181053 0.011078 16.344 < 2e-16 ***
## view3 0.274696 0.015238 18.027 < 2e-16 ***
## view4 0.337649 0.023231 14.534 < 2e-16 ***
## condition2 0.055359 0.063529 0.871 0.383548
## condition3 0.114851 0.058301 1.970 0.048861 *
## condition4 0.199371 0.058364 3.416 0.000637 ***
## condition5 0.279342 0.058767 4.753 2.02e-06 ***
## gradeAceptable -0.722114 0.279092 -2.587 0.009681 **
## gradeBuena -0.506090 0.279220 -1.813 0.069927 .
## gradeExcelente -0.071329 0.279833 -0.255 0.798804
## yr_renovated1 0.158466 0.011504 13.775 < 2e-16 ***
## lat 1.558626 0.016531 94.288 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2788 on 15014 degrees of freedom
## Multiple R-squared: 0.7201, Adjusted R-squared: 0.7198
## F-statistic: 2415 on 16 and 15014 DF, p-value: < 2.2e-16model3= lm(price~bathrooms+sqft_living+waterfront+view+condition+
grade+yr_built+lat, data=d_training) #modelo1 sin bedroom sin bedroom ni yr_renovated
summary(model3)#coeficientes que no son lógicos y R² = 0.7286 ##
## Call:
## lm(formula = price ~ bathrooms + sqft_living + waterfront + view +
## condition + grade + yr_built + lat, data = d_training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.20406 -0.17547 0.00073 0.16582 1.22462
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.500e+01 9.066e-01 -60.664 < 2e-16 ***
## bathrooms 5.774e-02 4.211e-03 13.710 < 2e-16 ***
## sqft_living 5.281e-01 8.657e-03 61.009 < 2e-16 ***
## waterfront1 3.771e-01 3.250e-02 11.604 < 2e-16 ***
## view1 2.075e-01 1.791e-02 11.585 < 2e-16 ***
## view2 1.528e-01 1.097e-02 13.922 < 2e-16 ***
## view3 2.440e-01 1.506e-02 16.198 < 2e-16 ***
## view4 3.189e-01 2.287e-02 13.943 < 2e-16 ***
## condition2 8.701e-02 6.254e-02 1.391 0.164156
## condition3 1.973e-01 5.746e-02 3.433 0.000598 ***
## condition4 2.380e-01 5.746e-02 4.141 3.48e-05 ***
## condition5 2.877e-01 5.784e-02 4.974 6.61e-07 ***
## gradeAceptable -6.492e-01 2.747e-01 -2.363 0.018132 *
## gradeBuena -4.032e-01 2.748e-01 -1.467 0.142422
## gradeExcelente 2.096e-02 2.754e-01 0.076 0.939333
## yr_built -2.713e-03 1.039e-04 -26.103 < 2e-16 ***
## lat 1.462e+00 1.673e-02 87.393 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2744 on 15014 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7286
## F-statistic: 2523 on 16 and 15014 DF, p-value: < 2.2e-16Teniendo en cuenta el R² y los coeficientes de las variables seleccionadas en cada modelo anterior nos quedaremos con el modelo denominado “model1a” :
model1a=lm(price~sqft_living+waterfront+view+condition+
+grade+yr_renovated+yr_built+lat, data=d_training) y veamos cómo de bien o mal precide el precio de las casas del conjunto train.
Colinealidad en el modelo Para ver si existe colinealidad entre las variables seleccionadas en nuestro modelo seleccionado vamos a usar el Factor de Inflación de la Varianza (VIF),
Los límites de referencia que se suelen emplear son:
- VIF = 1: Ausencia total de colinealidad.
- 1 < VIF < 5: La regresión puede verse afectada por cierta colinealidad.
- 5 < VIF < 10: Causa de preocupación.
vif(model1a)## GVIF Df GVIF^(1/(2*Df))
## sqft_living 1.785524 1 1.336235
## waterfront 1.506865 1 1.227544
## view 1.678938 4 1.066914
## condition 1.282599 4 1.031600
## grade 1.937298 3 1.116519
## yr_renovated 1.123988 1 1.060183
## yr_built 1.734619 1 1.317050
## lat 1.075865 1 1.037239Observamos que los valores se encuentran entre 1 y 2 y por consiguiente no nos indican la presencia de colinealidad entre variables.
Análisis de residuos
Una vez estimado el modelo de regresión y obtenido los residuos hay que comprobar si las hipótesis que se han utilizado para construirlo se pueden asumir como ciertas o no. Si no lo son, habrá que modificar el modelo para adaptarlo a los datos observados.
Para hacer el modelo se asumen 4 hipótesis: * 1 Linealidad: La relación entre la variable respuesta y las variables explicativas es lineal. * 2 Homocedasticidad: La variabilidad del error es constante. Esto es, el error sigue una distribución normal con varianza desconocida, pero todas iguales y constantes. * 3 Normalidad: Las perturbaciones (el error) siguen una distribución normal. * 4 Independencia: Las perturbaciones son independientes entre sí.
par(mfrow=c(2,2))
plot(model1a)
Linealidad
Se está analizando si existe una relación lineal entre la variable respuesta y las variables explicativas. Esto lo comprobamos en la primera de las 4 gráficas que tenemos.
Si se verifica la hipótesis de linealidad, esta gráfica debería de presentar una simetría respecto al eje horizontal. Como vemos en nuestro gráfico si se verifica la linealidad.
Normalidad
Para comprobar la normalidad de los residuos hay que hacer un gráfico cuantil-cuantil (un Q-Q plot) de los residuos, que se corresponde con la gŕafica de la esquina superior derecha.
A simple vista podŕiamos decir que a pesar de las colas si siguen una distribución normal. Comprobémoslo usando un contraste de normalidad:
residuos=resid(model1a)
ks.test(residuos, 'pnorm')##
## One-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.27828, p-value < 2.2e-16
## alternative hypothesis: two-sidedSe contrasta la hipótesis nula de que los residuos del modelo se distribuyen según una Normal. Al obtener un p.valor menor que 0.05 rechazamos H0 y por consiguiente se tiene que los residuos no siguen una distribución Normal. Sin embargo, debido a que los modelos lineales en general son robustos a la ligera falta de normalidad asumiremos la hipótesis de normalidad como cierta y proseguiremos con el estudio.
Homocedasticidad
Para comprobar la homocedasticidad hay que hacer un gráfico de dispersión cuantil-cuantil de los residuos. Lo ideal es que la recta sea parecida a y=0. Debemos obtener una línea recta horizontal. Si no fuese así lo que nos está mostrando es cierta dependencia entre la magnitud de los errores y la predicción.
El gŕafico inferior derecho es un gráfico de leverage, de apalancamiento, que sirve para comprobar si hay alguna observación que sea demasiado influyente en la construcción de los coeficientes del modelo. Si los puntos están agrupados y sin sobrepasar las curvas de nivel de la distancia de cook no habrá ningún problema. Si hay algún dato muy influyente, el gráfico te da un número indicando cuál es exactamente y tendríamos que estudiarlas para ver si son un error o un dato verdaderamente distinto al resto (un dato que sigue un patrón muy distinto al resto) e investigar su causa. El estudio de estos valores atípicos puede ser interesante porque pueden modificar los coeficientes del modelo
7. Predicción sobre los datos de testing. Evaluación del modelo
Una vez construido el modelo predictor, comprobemos como de bueno o malo o es. Esto lo haremos sobre el conjunto test que tenemos. Vamos a comprobar si con el modelo elegido obtenemos una buena predicción del precio.
Primero, realizaremos las mismas transformaciones, agrupaciones y categorizaciones que realizamos con los datos de train.
d_testing$price <- log(d_testing$price)
d_testing$bathrooms=ceiling(d_testing$bathrooms)
d_testing$sqft_living<- log(d_testing$sqft_living)
d_testing$grade <- cut(d_testing$grade, breaks = c(0,3,7,10,13), labels = c("Inacabada","Aceptable","Buena","Excelente"))
d_testing$basement=ifelse(d_testing$sqft_living-log(d_testing$sqft_above)==0,"0","1")
d_testing$basement=as.factor(d_testing$basement)
d_testing$yr_renovated=ifelse(d_testing$yr_renovated==0,"0", "1")
d_testing$yr_renovated=as.factor(d_testing$yr_renovated)Ahora veamos cómo predice el modelo:
pred=predict(model1a,d_testing[,-3])
New=as.data.frame(cbind(Prediccion=exp(pred),Observ=exp(d_testing[,3]), Dif=exp(pred)-exp(d_testing[,3])))
summary(New)## Prediccion Observ Dif
## Min. : 147591 Min. : 78000 Min. :-3443298
## 1st Qu.: 337479 1st Qu.: 321625 1st Qu.: -80093
## Median : 450798 Median : 450000 Median : -1071
## Mean : 515613 Mean : 538148 Mean : -22416
## 3rd Qu.: 629282 3rd Qu.: 645000 3rd Qu.: 71015
## Max. :4699766 Max. :7700000 Max. : 1514141
## NA's :22 NA's :22Obtenemos que el valor absoluto de la media de la diferencia entre el precio predicho y el observado es de 22416. Debido a que tenemos un error pequeño de predicción podemos considerar que tenemos un buen modelo.
plot(exp(d_testing$price),exp(pred),col="blue", xlab = "Precio según el modelo",
ylab = "Precio observado", main="Representación del precio predicho frente al observado")